home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / lapack / zhseqr.f < prev    next >
Text File  |  1996-07-19  |  15KB  |  462 lines
  1.       SUBROUTINE ZHSEQR( JOB, COMPZ, N, ILO, IHI, H, LDH, W, Z, LDZ,
  2.      $                   WORK, LWORK, INFO )
  3. *
  4. *  -- LAPACK routine (version 2.0) --
  5. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  6. *     Courant Institute, Argonne National Lab, and Rice University
  7. *     September 30, 1994
  8. *
  9. *     .. Scalar Arguments ..
  10.       CHARACTER          COMPZ, JOB
  11.       INTEGER            IHI, ILO, INFO, LDH, LDZ, LWORK, N
  12. *     ..
  13. *     .. Array Arguments ..
  14.       COMPLEX*16         H( LDH, * ), W( * ), WORK( * ), Z( LDZ, * )
  15. *     ..
  16. *
  17. *  Purpose
  18. *  =======
  19. *
  20. *  ZHSEQR computes the eigenvalues of a complex upper Hessenberg
  21. *  matrix H, and, optionally, the matrices T and Z from the Schur
  22. *  decomposition H = Z T Z**H, where T is an upper triangular matrix
  23. *  (the Schur form), and Z is the unitary matrix of Schur vectors.
  24. *
  25. *  Optionally Z may be postmultiplied into an input unitary matrix Q,
  26. *  so that this routine can give the Schur factorization of a matrix A
  27. *  which has been reduced to the Hessenberg form H by the unitary
  28. *  matrix Q:  A = Q*H*Q**H = (QZ)*T*(QZ)**H.
  29. *
  30. *  Arguments
  31. *  =========
  32. *
  33. *  JOB     (input) CHARACTER*1
  34. *          = 'E': compute eigenvalues only;
  35. *          = 'S': compute eigenvalues and the Schur form T.
  36. *
  37. *  COMPZ   (input) CHARACTER*1
  38. *          = 'N': no Schur vectors are computed;
  39. *          = 'I': Z is initialized to the unit matrix and the matrix Z
  40. *                 of Schur vectors of H is returned;
  41. *          = 'V': Z must contain an unitary matrix Q on entry, and
  42. *                 the product Q*Z is returned.
  43. *
  44. *  N       (input) INTEGER
  45. *          The order of the matrix H.  N >= 0.
  46. *
  47. *  ILO     (input) INTEGER
  48. *  IHI     (input) INTEGER
  49. *          It is assumed that H is already upper triangular in rows
  50. *          and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally
  51. *          set by a previous call to ZGEBAL, and then passed to CGEHRD
  52. *          when the matrix output by ZGEBAL is reduced to Hessenberg
  53. *          form. Otherwise ILO and IHI should be set to 1 and N
  54. *          respectively.
  55. *          1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
  56. *
  57. *  H       (input/output) COMPLEX*16 array, dimension (LDH,N)
  58. *          On entry, the upper Hessenberg matrix H.
  59. *          On exit, if JOB = 'S', H contains the upper triangular matrix
  60. *          T from the Schur decomposition (the Schur form). If
  61. *          JOB = 'E', the contents of H are unspecified on exit.
  62. *
  63. *  LDH     (input) INTEGER
  64. *          The leading dimension of the array H. LDH >= max(1,N).
  65. *
  66. *  W       (output) COMPLEX*16 array, dimension (N)
  67. *          The computed eigenvalues. If JOB = 'S', the eigenvalues are
  68. *          stored in the same order as on the diagonal of the Schur form
  69. *          returned in H, with W(i) = H(i,i).
  70. *
  71. *  Z       (input/output) COMPLEX*16 array, dimension (LDZ,N)
  72. *          If COMPZ = 'N': Z is not referenced.
  73. *          If COMPZ = 'I': on entry, Z need not be set, and on exit, Z
  74. *          contains the unitary matrix Z of the Schur vectors of H.
  75. *          If COMPZ = 'V': on entry Z must contain an N-by-N matrix Q,
  76. *          which is assumed to be equal to the unit matrix except for
  77. *          the submatrix Z(ILO:IHI,ILO:IHI); on exit Z contains Q*Z.
  78. *          Normally Q is the unitary matrix generated by ZUNGHR after
  79. *          the call to ZGEHRD which formed the Hessenberg matrix H.
  80. *
  81. *  LDZ     (input) INTEGER
  82. *          The leading dimension of the array Z.
  83. *          LDZ >= max(1,N) if COMPZ = 'I' or 'V'; LDZ >= 1 otherwise.
  84. *
  85. *  WORK    (workspace) COMPLEX*16 array, dimension (N)
  86. *
  87. *  LWORK   (input) INTEGER
  88. *          This argument is currently redundant.
  89. *
  90. *  INFO    (output) INTEGER
  91. *          = 0:  successful exit
  92. *          < 0:  if INFO = -i, the i-th argument had an illegal value
  93. *          > 0:  if INFO = i, ZHSEQR failed to compute all the
  94. *                eigenvalues in a total of 30*(IHI-ILO+1) iterations;
  95. *                elements 1:ilo-1 and i+1:n of W contain those
  96. *                eigenvalues which have been successfully computed.
  97. *
  98. *  =====================================================================
  99. *
  100. *     .. Parameters ..
  101.       COMPLEX*16         ZERO, ONE
  102.       PARAMETER          ( ZERO = ( 0.0D+0, 0.0D+0 ),
  103.      $                   ONE = ( 1.0D+0, 0.0D+0 ) )
  104.       DOUBLE PRECISION   RZERO, RONE, CONST
  105.       PARAMETER          ( RZERO = 0.0D+0, RONE = 1.0D+0,
  106.      $                   CONST = 1.5D+0 )
  107.       INTEGER            NSMAX, LDS
  108.       PARAMETER          ( NSMAX = 15, LDS = NSMAX )
  109. *     ..
  110. *     .. Local Scalars ..
  111.       LOGICAL            INITZ, WANTT, WANTZ
  112.       INTEGER            I, I1, I2, IERR, II, ITEMP, ITN, ITS, J, K, L,
  113.      $                   MAXB, NH, NR, NS, NV
  114.       DOUBLE PRECISION   OVFL, RTEMP, SMLNUM, TST1, ULP, UNFL
  115.       COMPLEX*16         CDUM, TAU, TEMP
  116. *     ..
  117. *     .. Local Arrays ..
  118.       DOUBLE PRECISION   RWORK( 1 )
  119.       COMPLEX*16         S( LDS, NSMAX ), V( NSMAX+1 ), VV( NSMAX+1 )
  120. *     ..
  121. *     .. External Functions ..
  122.       LOGICAL            LSAME
  123.       INTEGER            ILAENV, IZAMAX
  124.       DOUBLE PRECISION   DLAMCH, DLAPY2, ZLANHS
  125.       EXTERNAL           LSAME, ILAENV, IZAMAX, DLAMCH, DLAPY2, ZLANHS
  126. *     ..
  127. *     .. External Subroutines ..
  128.       EXTERNAL           DLABAD, XERBLA, ZCOPY, ZDSCAL, ZGEMV, ZLACPY,
  129.      $                   ZLAHQR, ZLARFG, ZLARFX, ZLASET, ZSCAL
  130. *     ..
  131. *     .. Intrinsic Functions ..
  132.       INTRINSIC          ABS, DBLE, DCONJG, DIMAG, MAX, MIN
  133. *     ..
  134. *     .. Statement Functions ..
  135.       DOUBLE PRECISION   CABS1
  136. *     ..
  137. *     .. Statement Function definitions ..
  138.       CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
  139. *     ..
  140. *     .. Executable Statements ..
  141. *
  142. *     Decode and test the input parameters
  143. *
  144.       WANTT = LSAME( JOB, 'S' )
  145.       INITZ = LSAME( COMPZ, 'I' )
  146.       WANTZ = INITZ .OR. LSAME( COMPZ, 'V' )
  147. *
  148.       INFO = 0
  149.       IF( .NOT.LSAME( JOB, 'E' ) .AND. .NOT.WANTT ) THEN
  150.          INFO = -1
  151.       ELSE IF( .NOT.LSAME( COMPZ, 'N' ) .AND. .NOT.WANTZ ) THEN
  152.          INFO = -2
  153.       ELSE IF( N.LT.0 ) THEN
  154.          INFO = -3
  155.       ELSE IF( ILO.LT.1 .OR. ILO.GT.MAX( 1, N ) ) THEN
  156.          INFO = -4
  157.       ELSE IF( IHI.LT.MIN( ILO, N ) .OR. IHI.GT.N ) THEN
  158.          INFO = -5
  159.       ELSE IF( LDH.LT.MAX( 1, N ) ) THEN
  160.          INFO = -7
  161.       ELSE IF( LDZ.LT.1 .OR. WANTZ .AND. LDZ.LT.MAX( 1, N ) ) THEN
  162.          INFO = -10
  163.       END IF
  164.       IF( INFO.NE.0 ) THEN
  165.          CALL XERBLA( 'ZHSEQR', -INFO )
  166.          RETURN
  167.       END IF
  168. *
  169. *     Initialize Z, if necessary
  170. *
  171.       IF( INITZ )
  172.      $   CALL ZLASET( 'Full', N, N, ZERO, ONE, Z, LDZ )
  173. *
  174. *     Store the eigenvalues isolated by ZGEBAL.
  175. *
  176.       DO 10 I = 1, ILO - 1
  177.          W( I ) = H( I, I )
  178.    10 CONTINUE
  179.       DO 20 I = IHI + 1, N
  180.          W( I ) = H( I, I )
  181.    20 CONTINUE
  182. *
  183. *     Quick return if possible.
  184. *
  185.       IF( N.EQ.0 )
  186.      $   RETURN
  187.       IF( ILO.EQ.IHI ) THEN
  188.          W( ILO ) = H( ILO, ILO )
  189.          RETURN
  190.       END IF
  191. *
  192. *     Set rows and columns ILO to IHI to zero below the first
  193. *     subdiagonal.
  194. *
  195.       DO 40 J = ILO, IHI - 2
  196.          DO 30 I = J + 2, N
  197.             H( I, J ) = ZERO
  198.    30    CONTINUE
  199.    40 CONTINUE
  200.       NH = IHI - ILO + 1
  201. *
  202. *     I1 and I2 are the indices of the first row and last column of H
  203. *     to which transformations must be applied. If eigenvalues only are
  204. *     being computed, I1 and I2 are re-set inside the main loop.
  205. *
  206.       IF( WANTT ) THEN
  207.          I1 = 1
  208.          I2 = N
  209.       ELSE
  210.          I1 = ILO
  211.          I2 = IHI
  212.       END IF
  213. *
  214. *     Ensure that the subdiagonal elements are real.
  215. *
  216.       DO 50 I = ILO + 1, IHI
  217.          TEMP = H( I, I-1 )
  218.          IF( DIMAG( TEMP ).NE.RZERO ) THEN
  219.             RTEMP = DLAPY2( DBLE( TEMP ), DIMAG( TEMP ) )
  220.             H( I, I-1 ) = RTEMP
  221.             TEMP = TEMP / RTEMP
  222.             IF( I2.GT.I )
  223.      $         CALL ZSCAL( I2-I, DCONJG( TEMP ), H( I, I+1 ), LDH )
  224.             CALL ZSCAL( I-I1, TEMP, H( I1, I ), 1 )
  225.             IF( I.LT.IHI )
  226.      $         H( I+1, I ) = TEMP*H( I+1, I )
  227.             IF( WANTZ )
  228.      $         CALL ZSCAL( NH, TEMP, Z( ILO, I ), 1 )
  229.          END IF
  230.    50 CONTINUE
  231. *
  232. *     Determine the order of the multi-shift QR algorithm to be used.
  233. *
  234.       NS = ILAENV( 4, 'ZHSEQR', JOB // COMPZ, N, ILO, IHI, -1 )
  235.       MAXB = ILAENV( 8, 'ZHSEQR', JOB // COMPZ, N, ILO, IHI, -1 )
  236.       IF( NS.LE.1 .OR. NS.GT.NH .OR. MAXB.GE.NH ) THEN
  237. *
  238. *        Use the standard double-shift algorithm
  239. *
  240.          CALL ZLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILO, IHI, Z,
  241.      $                LDZ, INFO )
  242.          RETURN
  243.       END IF
  244.       MAXB = MAX( 2, MAXB )
  245.       NS = MIN( NS, MAXB, NSMAX )
  246. *
  247. *     Now 1 < NS <= MAXB < NH.
  248. *
  249. *     Set machine-dependent constants for the stopping criterion.
  250. *     If norm(H) <= sqrt(OVFL), overflow should not occur.
  251. *
  252.       UNFL = DLAMCH( 'Safe minimum' )
  253.       OVFL = RONE / UNFL
  254.       CALL DLABAD( UNFL, OVFL )
  255.       ULP = DLAMCH( 'Precision' )
  256.       SMLNUM = UNFL*( NH / ULP )
  257. *
  258. *     ITN is the total number of multiple-shift QR iterations allowed.
  259. *
  260.       ITN = 30*NH
  261. *
  262. *     The main loop begins here. I is the loop index and decreases from
  263. *     IHI to ILO in steps of at most MAXB. Each iteration of the loop
  264. *     works with the active submatrix in rows and columns L to I.
  265. *     Eigenvalues I+1 to IHI have already converged. Either L = ILO, or
  266. *     H(L,L-1) is negligible so that the matrix splits.
  267. *
  268.       I = IHI
  269.    60 CONTINUE
  270.       IF( I.LT.ILO )
  271.      $   GO TO 180
  272. *
  273. *     Perform multiple-shift QR iterations on rows and columns ILO to I
  274. *     until a submatrix of order at most MAXB splits off at the bottom
  275. *     because a subdiagonal element has become negligible.
  276. *
  277.       L = ILO
  278.       DO 160 ITS = 0, ITN
  279. *
  280. *        Look for a single small subdiagonal element.
  281. *
  282.          DO 70 K = I, L + 1, -1
  283.             TST1 = CABS1( H( K-1, K-1 ) ) + CABS1( H( K, K ) )
  284.             IF( TST1.EQ.RZERO )
  285.      $         TST1 = ZLANHS( '1', I-L+1, H( L, L ), LDH, RWORK )
  286.             IF( ABS( DBLE( H( K, K-1 ) ) ).LE.MAX( ULP*TST1, SMLNUM ) )
  287.      $         GO TO 80
  288.    70    CONTINUE
  289.    80    CONTINUE
  290.          L = K
  291.          IF( L.GT.ILO ) THEN
  292. *
  293. *           H(L,L-1) is negligible.
  294. *
  295.             H( L, L-1 ) = ZERO
  296.          END IF
  297. *
  298. *        Exit from loop if a submatrix of order <= MAXB has split off.
  299. *
  300.          IF( L.GE.I-MAXB+1 )
  301.      $      GO TO 170
  302. *
  303. *        Now the active submatrix is in rows and columns L to I. If
  304. *        eigenvalues only are being computed, only the active submatrix
  305. *        need be transformed.
  306. *
  307.          IF( .NOT.WANTT ) THEN
  308.             I1 = L
  309.             I2 = I
  310.          END IF
  311. *
  312.          IF( ITS.EQ.20 .OR. ITS.EQ.30 ) THEN
  313. *
  314. *           Exceptional shifts.
  315. *
  316.             DO 90 II = I - NS + 1, I
  317.                W( II ) = CONST*( ABS( DBLE( H( II, II-1 ) ) )+
  318.      $                   ABS( DBLE( H( II, II ) ) ) )
  319.    90       CONTINUE
  320.          ELSE
  321. *
  322. *           Use eigenvalues of trailing submatrix of order NS as shifts.
  323. *
  324.             CALL ZLACPY( 'Full', NS, NS, H( I-NS+1, I-NS+1 ), LDH, S,
  325.      $                   LDS )
  326.             CALL ZLAHQR( .FALSE., .FALSE., NS, 1, NS, S, LDS,
  327.      $                   W( I-NS+1 ), 1, NS, Z, LDZ, IERR )
  328.             IF( IERR.GT.0 ) THEN
  329. *
  330. *              If ZLAHQR failed to compute all NS eigenvalues, use the
  331. *              unconverged diagonal elements as the remaining shifts.
  332. *
  333.                DO 100 II = 1, IERR
  334.                   W( I-NS+II ) = S( II, II )
  335.   100          CONTINUE
  336.             END IF
  337.          END IF
  338. *
  339. *        Form the first column of (G-w(1)) (G-w(2)) . . . (G-w(ns))
  340. *        where G is the Hessenberg submatrix H(L:I,L:I) and w is
  341. *        the vector of shifts (stored in W). The result is
  342. *        stored in the local array V.
  343. *
  344.          V( 1 ) = ONE
  345.          DO 110 II = 2, NS + 1
  346.             V( II ) = ZERO
  347.   110    CONTINUE
  348.          NV = 1
  349.          DO 130 J = I - NS + 1, I
  350.             CALL ZCOPY( NV+1, V, 1, VV, 1 )
  351.             CALL ZGEMV( 'No transpose', NV+1, NV, ONE, H( L, L ), LDH,
  352.      $                  VV, 1, -W( J ), V, 1 )
  353.             NV = NV + 1
  354. *
  355. *           Scale V(1:NV) so that max(abs(V(i))) = 1. If V is zero,
  356. *           reset it to the unit vector.
  357. *
  358.             ITEMP = IZAMAX( NV, V, 1 )
  359.             RTEMP = CABS1( V( ITEMP ) )
  360.             IF( RTEMP.EQ.RZERO ) THEN
  361.                V( 1 ) = ONE
  362.                DO 120 II = 2, NV
  363.                   V( II ) = ZERO
  364.   120          CONTINUE
  365.             ELSE
  366.                RTEMP = MAX( RTEMP, SMLNUM )
  367.                CALL ZDSCAL( NV, RONE / RTEMP, V, 1 )
  368.             END IF
  369.   130    CONTINUE
  370. *
  371. *        Multiple-shift QR step
  372. *
  373.          DO 150 K = L, I - 1
  374. *
  375. *           The first iteration of this loop determines a reflection G
  376. *           from the vector V and applies it from left and right to H,
  377. *           thus creating a nonzero bulge below the subdiagonal.
  378. *
  379. *           Each subsequent iteration determines a reflection G to
  380. *           restore the Hessenberg form in the (K-1)th column, and thus
  381. *           chases the bulge one step toward the bottom of the active
  382. *           submatrix. NR is the order of G.
  383. *
  384.             NR = MIN( NS+1, I-K+1 )
  385.             IF( K.GT.L )
  386.      $         CALL ZCOPY( NR, H( K, K-1 ), 1, V, 1 )
  387.             CALL ZLARFG( NR, V( 1 ), V( 2 ), 1, TAU )
  388.             IF( K.GT.L ) THEN
  389.                H( K, K-1 ) = V( 1 )
  390.                DO 140 II = K + 1, I
  391.                   H( II, K-1 ) = ZERO
  392.   140          CONTINUE
  393.             END IF
  394.             V( 1 ) = ONE
  395. *
  396. *           Apply G' from the left to transform the rows of the matrix
  397. *           in columns K to I2.
  398. *
  399.             CALL ZLARFX( 'Left', NR, I2-K+1, V, DCONJG( TAU ),
  400.      $                   H( K, K ), LDH, WORK )
  401. *
  402. *           Apply G from the right to transform the columns of the
  403. *           matrix in rows I1 to min(K+NR,I).
  404. *
  405.             CALL ZLARFX( 'Right', MIN( K+NR, I )-I1+1, NR, V, TAU,
  406.      $                   H( I1, K ), LDH, WORK )
  407. *
  408.             IF( WANTZ ) THEN
  409. *
  410. *              Accumulate transformations in the matrix Z
  411. *
  412.                CALL ZLARFX( 'Right', NH, NR, V, TAU, Z( ILO, K ), LDZ,
  413.      $                      WORK )
  414.             END IF
  415.   150    CONTINUE
  416. *
  417. *        Ensure that H(I,I-1) is real.
  418. *
  419.          TEMP = H( I, I-1 )
  420.          IF( DIMAG( TEMP ).NE.RZERO ) THEN
  421.             RTEMP = DLAPY2( DBLE( TEMP ), DIMAG( TEMP ) )
  422.             H( I, I-1 ) = RTEMP
  423.             TEMP = TEMP / RTEMP
  424.             IF( I2.GT.I )
  425.      $         CALL ZSCAL( I2-I, DCONJG( TEMP ), H( I, I+1 ), LDH )
  426.             CALL ZSCAL( I-I1, TEMP, H( I1, I ), 1 )
  427.             IF( WANTZ ) THEN
  428.                CALL ZSCAL( NH, TEMP, Z( ILO, I ), 1 )
  429.             END IF
  430.          END IF
  431. *
  432.   160 CONTINUE
  433. *
  434. *     Failure to converge in remaining number of iterations
  435. *
  436.       INFO = I
  437.       RETURN
  438. *
  439.   170 CONTINUE
  440. *
  441. *     A submatrix of order <= MAXB in rows and columns L to I has split
  442. *     off. Use the double-shift QR algorithm to handle it.
  443. *
  444.       CALL ZLAHQR( WANTT, WANTZ, N, L, I, H, LDH, W, ILO, IHI, Z, LDZ,
  445.      $             INFO )
  446.       IF( INFO.GT.0 )
  447.      $   RETURN
  448. *
  449. *     Decrement number of remaining iterations, and return to start of
  450. *     the main loop with a new value of I.
  451. *
  452.       ITN = ITN - ITS
  453.       I = L - 1
  454.       GO TO 60
  455. *
  456.   180 CONTINUE
  457.       RETURN
  458. *
  459. *     End of ZHSEQR
  460. *
  461.       END
  462.